import numpy as np
from scipy import integrate
from scipy.optimize import brentq,fsolve, root,bisect
from sympy import diff, log
import sympy as sp
import matplotlib.pyplot as plt
import pandas as pd
from . import PE_cal 
import traceback

# 相互作用系数 La Ce Gd Y Nd
sys = ['Mg_Al','Mg_La','Mg_Ce','Mg_Gd','Mg_Y','Mg_Nd','Al_La','Al_Ce','Al_Gd','Al_Y','Al_Nd']

def parameters(system,T):
    try:
            
        L12,theta,M,rou=[],[],[],[]
        # theta M rou
        Mg = [(577-0.26*(T-923))/1000, 0.024305, 1740]
        Al = [(871-0.155*(T-933))/1000, 0.026982, 2700]
        Zn = [(789-0.21*(T-693))/1000, 0.06539, 7140]
        La = [(728-0.1*(T-1193))/1000, 0.138906, 6150]
        Ce = [(707-0.078*(T-1077))/1000, 0.140116, 6770]
        Gd = [(664-0.058*(T-1625))/1000, 0.15725, 7900]
        Y = [(872-0.086*(T-1793))/1000, 0.088906, 4470]
        Nd = [(687-0.0885*(T-1297))/1000, 0.14424, 7010]
        # 确保T是浮点数
        T = float(T)
        
        if system == 'Mg_Al':
            L12 = [-12000+T*8.566, 1890-3*T, 2000, 0]
            theta = [Mg[0], Al[0]]
            M = [Mg[1], Al[1]]
            rou = [Mg[2], Al[2]]
        elif system == 'Mg_La':
            L12 = [-33288.5+T*8.9095, -(41081-29.4287*T), -12051.9, 0]
            theta = [Mg[0], La[0]]
            M = [Mg[1], La[1]]
            rou = [Mg[2], La[2]]
        elif system == 'Mg_Ce':
            L12 = [-39381.2+T*16.3405, -25338.6, -15106.9, 0]
            theta = [Mg[0], Ce[0]]
            M = [Mg[1], Ce[1]]
            rou = [Mg[2], Ce[2]]
        elif system == 'Mg_Gd':
            L12 = [-36681.3+T*16.2484, -(34233.8-10.7783*T), -7352.9, 0]
            theta = [Mg[0], Gd[0]]
            M = [Mg[1], Gd[1]]
            rou = [Mg[2], Gd[2]]
        elif system == 'Mg_Y':
            L12 = [-40917.97+T*22.83, -18685.23+10.24*T, 1076.67-6.05*T, 0]
            theta = [Mg[0], Y[0]]
            M = [Mg[1], Y[1]]
            rou = [Mg[2], Y[2]]
        elif system == 'Mg_Nd':
            L12 = [-34308.7+T*10.9, 14565.6, -16479, 0]
            theta = [Mg[0], Nd[0]]
            M = [Mg[1], Nd[1]]
            rou = [Mg[2], Nd[2]]
        elif system == 'Al_La':
            L12 = [-89252, 0, 0, 0]
            theta = [Al[0], La[0]]
            M = [Al[1], La[1]]
            rou = [Al[2], La[2]]
        elif system == 'Al_Ce':
            L12 = [-141000+43.82*T, -55600+23.68*T, -1686, 0]
            theta = [Al[0], Ce[0]]
            M = [Al[1], Ce[1]]
            rou = [Al[2], Ce[2]]
        elif system == 'Al_Gd':
            L12 = [-166500+52.36*T, -23790, 18520, 0]
            theta = [Al[0], Gd[0]]
            M = [Al[1], Gd[1]]
            rou = [Al[2], Gd[2]]
        elif system == 'Al_Y':
            L12 = [-160876.36+5*T, -32000+7.56*T, 32000-6.53*T, 0]
            theta = [Al[0], Y[0]]
            M = [Al[1], Y[1]]
            rou = [Al[2], Y[2]]
        elif system == 'Al_Nd':
            L12 = [-145100+44.97*T, -28030, -6708, 0]
            theta = [Al[0], Nd[0]]
            M = [Al[1], Nd[1]]
            rou = [Al[2], Nd[2]]
        elif system == 'Mg_Zn':
            L12 = [-77729.2+T*680.523-95*T*log(T)+0.04*T**2,1890-3*T,2000,0]
            theta = [Mg[0],Zn[0]]
            M = [Mg[1],Zn[1]]
            rou = [Mg[2],Zn[2]]
        
        elif system == 'Al_Zn':
            L12 = [10465.55-3.39259*T,0,0,0]
            theta = [Al[0],Zn[0]]
            M = [Al[1],Zn[1]]
            rou = [Al[2],Zn[2]]

        elif system == 'Zn_La':
            L12 = [-89252,0,0,0]
            theta = [Zn[0],La[0]]
            M = [Zn[1],La[1]]
            rou = [Zn[2],La[2]]
        elif system == 'Zn_Ce':
            L12 = [-141000+43.82*T,-55600+23.68*T,-1686,0]
            theta = [Zn[0],Ce[0]]
            M = [Zn[1],Ce[1]]
            rou = [Zn[2],Ce[2]]
        elif system == 'Zn_Gd':
            L12 = [-166500+52.36*T,-23790,18520,0]
            theta = [Zn[0],Gd[0]]
            M = [Zn[1],Gd[1]]
            rou = [Zn[2],Gd[2]]
        elif system == 'Zn_Y':
            L12 = [-160876.36+5*T,-32000+7.56*T,32000-6.53*T,0]
            theta = [Zn[0],Y[0]]
            M = [Zn[1],Y[1]]
            rou = [Zn[2],Y[2]]
        elif system == 'Zn_Nd':
            L12 = [-145100+44.97*T,-28030,-6708,0]
            theta = [Zn[0],Nd[0]]
            M = [Zn[1],Nd[1]]
            rou = [Zn[2],Nd[2]]
        else:
            raise ValueError(f"Unsupported system: {system}")

        # 确保所有参数都是浮点数
        L12 = [float(x) for x in L12]
        theta = [float(x) for x in theta]
        M = [float(x) for x in M]
        rou = [float(x) for x in rou]

        return L12, theta, M, rou
        
    except Exception as e:
        print(f"Error in parameters function: {str(e)}")
        print(f"Stack trace: {traceback.format_exc()}")
        raise

# 超额黏度
def parameter_vicosity(system,T):

    #Mg-Al-Si-Fe-Cu
    #Mg-Al-La

    L12,V = [],[]

    #纯组元黏度
    Mg=0.04809*2.7182818**(3003.420496/T)
    Al=0.15957*2.7182818**(2000.740026/T)
    Si=0.19374*2.7182818**(1823.187981/T)
    Fe=0.19014*2.7182818**(6205.43843/T)
    Cu=0.42817*2.7182818**(2998.556651/T)
    La=1.297*2.7182818**(858.67/T)

    if system == 'Al_Cu':
        L12=[4.24870,-16.75004,2.86590,9.80701]
        V=[Al,Cu]
    if system == 'Al_Si':
        L12 = [-0.825619841, -1.79026909, 0.418292707,0]
        V=[Al,Si]
    if system == 'Si_Cu':
        L12 = [-3.27440891,-4.12012308,10.0780159,6.82682157]
        V = [Si,Cu]
    if system == 'Mg_Si':
        L12=[3.35157962,4.33610175,0,0]
        V=[Mg,Si]
    if system == 'Mg_Cu':
        L12 = [-4.59329456,2.76068945, -2.13354732,0]
        V=[Mg,Cu]
    if system == 'Mg_Al':
        L12 = [0.897678, 0.829784,0.544635,0]
        V=[Mg,Al]
    if system == 'Al_Fe':
        L12=[137-0.088*T,-197+0.12*T,-0.35,1.18]
        V=[Al,Fe]
    if system == 'Si_Fe':
        L12 = [-3.818E-06*T**2+1.771E-02*T-22.08,-(6.781E-06*T**2-3.070E-02*T+35.86), 3.947E-06*T**2-1.758E-02*T+20.01,-(-1.098E-06*T**2+5.029E-03*T-6.028)]
        V=[Si,Fe]
    if system == 'Fe_Cu':
        L12 = [-1.178E-05*T**2+5.243E-02*T-66.56, -(6.696E-06*T**2-3.031E-02*T+35.54),-4.021E-06*T**2+1.797E-02*T-2.064,-(2.017E-06*T**2-8.945E-03*T+10.1577)]
        V=[Fe,Cu]
    if system == 'Mg_La':
        L12 = [-0.675, 1.239,0.233,0]
        V=[Mg,La]
    if system == 'Al_La':
        L12 = [-1.314, 0.567,0.349,0]
        V=[Al,La]

    if system == 'Mg_Fe':
        L12 = [0.932,-4.786,1.038,1.127]
        V=[Mg,Fe]

    return L12,V


def butler(T,beta,l_12,theta,M,rou,x1):
    try:
        # 确保所有输入都是浮点数
        T = float(T)
        beta = float(beta)
        x1 = float(x1)
        
        # 确保列表中的元素都是浮点数，并检查长度
        l_12 = [float(x) for x in l_12]
        if len(l_12) < 4:
            l_12.extend([0.0] * (4 - len(l_12)))  # 确保l_12有4个元素
            
        theta = [float(x) for x in theta]
        M = [float(x) for x in M]
        rou = [float(x) for x in rou]
        
        # 打印调试信息
        print(f"Input parameters after conversion:")
        print(f"T: {T}, beta: {beta}, x1: {x1}")
        print(f"l_12: {l_12}")
        print(f"theta: {theta}")
        print(f"M: {M}")
        print(f"rou: {rou}")
        
        R = 8.314
        N0 = 6.02E+23
        # 二元 体系
        L = 1.091
        # 使用numpy的power函数替代pow
        A = [L * np.power(N0,1/3) * np.power(M[0]/rou[0], 2/3), 
             L * np.power(N0, 1/3) * np.power(M[1]/rou[1], 2/3)]
        
        # 打印中间计算结果
        print(f"Intermediate calculations:")
        print(f"A: {A}")
        
        # deltaG 表达式
        x2 = 1-x1
        # 使用numpy的power函数
        deltaG = x1 * x2 * (l_12[0] + l_12[1] * (x1 - x2) + 
                           l_12[2] * np.power(x1-x2, 2) + 
                           l_12[3] * np.power(x1-x2, 3))
        
        deltaG_B1 = deltaG + (1-x1)*((1-2*x1)*(l_12[0]+l_12[1]*(2*x1-1)+
                        l_12[2]*np.power(2*x1-1,2)+
                        l_12[3]*np.power(2*x1-1,3))+
                x1*(1-x1)*(2*l_12[1]+4*l_12[2]*(2*x1-1)+
                        6*l_12[3]*np.power(2*x1-1,2)))
        
        deltaG_B2 = deltaG + (1-x2)*((1-2*x2)*(l_12[0]+l_12[1]*(1-2*x2)+
                        l_12[2]*np.power(2*x2-1,2)+
                        l_12[3]*np.power(2*x2-1,3))+
                x2*(1-x2)*(-2*l_12[1]-4*l_12[2]*(1-2*x2)-
                        6*l_12[3]*np.power(1-2*x2,2)))

        # print(f"deltaG calculations:")
        # print(f"deltaG: {deltaG}")
        # print(f"deltaG_B1: {deltaG_B1}")
        # print(f"deltaG_B2: {deltaG_B2}")

        def equation(x):
            x = float(x)
            # 添加数值稳定性检查
            if x <= 0 or x >= 1:
                return float('inf')
                
            # 使用numpy的log函数，并添加数值稳定性检查
            log_x_x1 = np.log(np.maximum(x / x1, 1e-300))
            log_1x_x2 = np.log(np.maximum((1 - x) / x2, 1e-300))
            
            term1 = theta[0] + (R * T * log_x_x1 + (beta-1) * deltaG_B1) / A[0]
            term2 = theta[1] + (R * T * log_1x_x2 + (beta-1) * deltaG_B2) / A[1]
            result = term1 - term2
            
            print(f"Equation evaluation at x={x}: term1={term1}, term2={term2}, result={result}")
            return result

        # 使用更安全的边界值
        try:
            solution = bisect(equation, 1e-4, 1 - 1e-4)
            print(f"Solution found: {solution}")
        except ValueError as e:
            print(f"Error in bisect: {str(e)}")
            raise

        # 使用numpy的log函数，并添加数值稳定性检查
        log_solution_x1 = np.log(np.maximum(solution / x1, 1e-300))
        theta12 = theta[0] + (R * T * log_solution_x1 + (beta-1) * deltaG_B1) / A[0]

        # 限制输出的小数位数
        theta12 = round(theta12, 6)
        print(f"Final result: {theta12}")
        return theta12
        
    except Exception as e:
        print(f"Error in butler function: {str(e)}")
        print(f"Input parameters: T={T}, beta={beta}, l_12={l_12}, theta={theta}, M={M}, rou={rou}, x1={x1}")
        print(f"Stack trace: {traceback.format_exc()}")
        raise


def calculate_butler_3(T,beta,l_12,l_23,l_13,theta,M,rou,x1,x2,x3):

    l_31 = [l_13[0], -l_13[1], l_13[2], -l_13[3]]

    # Define sigma functions
    def sigma_E12(x):
        return x * (1 - x) * (l_12[0] + l_12[1] * (2 * x - 1) + l_12[2] * (2 * x - 1) ** 2 + l_12[3]*pow((2*x-1),3))
    def sigma_E13(x):
        return x * (1 - x) * (l_31[0] - l_31[1] * (2 * x - 1) + l_31[2] * (2 * x - 1) ** 2 - l_31[3]*pow((2*x-1),3))
    def sigma_E21(x):
        return x * (1 - x) * (l_12[0] - l_12[1] * (2 * x - 1) + l_12[2] * (2 * x - 1) ** 2 - l_12[3]*pow((2*x-1),3))
    def sigma_E23(x):
        return x * (1 - x) * (l_23[0] + l_23[1] * (2 * x - 1) + l_23[2] * (2 * x - 1) ** 2 + l_23[3]*pow((2*x-1),3))
    def sigma_E31(x):
        return x * (1 - x) * (l_31[0] + l_31[1] * (2 * x - 1) + l_31[2] * (2 * x - 1) ** 2 + l_31[3]*pow((2*x-1),3))
    def sigma_E32(x):
        return x * (1 - x) * (l_23[0] - l_23[1] * (2 * x - 1) + l_23[2] * (2 * x - 1) ** 2 - l_23[3]*pow((2*x-1),3))

    # Calculate eta values
    def integrand1(x):
        return (sigma_E12(x) - sigma_E13(x)) ** 2
    def integrand2(x):
        return (sigma_E21(x) - sigma_E23(x)) ** 2
    def integrand3(x):
        return (sigma_E31(x) - sigma_E32(x)) ** 2

    eta1 = integrate.quad(integrand1, 0, 1)[0]
    eta2 = integrate.quad(integrand2, 0, 1)[0]
    eta3 = integrate.quad(integrand3, 0, 1)[0]

    # Calculate epsilon values
    epsilon12 = eta1 / (eta1 + eta2)
    epsilon23 = eta2 / (eta2 + eta3)
    epsilon31 = eta3 / (eta3 + eta1)
    # print(eta1, eta2, eta3, epsilon12, epsilon23, epsilon31)

    # Calculate X values
    X1 = x1 + x3 * epsilon12
    X2 = x2 + x1 * epsilon23
    X3 = x3 + x2 * epsilon31

    # Calculate w values
    w12 = (x1 * x2) / (X1 * (1 - X1))
    w23 = (x2 * x3) / (X2 * (1 - X2))
    w31 = (x3 * x1) / (X3 * (1 - X3))

    # Calculate G
    G = w12 * sigma_E12(X1) + w23 * sigma_E23(X2) + w31 * sigma_E31(X3)

    # Calculate partial derivatives using finite differences
    h = 1e-10

    def calculate_G(x1_val, x2_val, x3_val):
        # Calculate new X values based on new x values
        X1_new = x1_val + x3_val * epsilon12
        X2_new = x2_val + x1_val * epsilon23
        X3_new = x3_val + x2_val * epsilon31

        w12_new = (x1_val * x2_val) / (X1_new * (1 - X1_new))
        w23_new = (x2_val * x3_val) / (X2_new * (1 - X2_new))
        w31_new = (x3_val * x1_val) / (X3_new * (1 - X3_new))

        return w12_new * sigma_E12(X1_new) + w23_new * sigma_E23(X2_new) + w31_new * sigma_E31(X3_new)

    # Calculate partial derivatives
    dG1 = (calculate_G(x1 + h, x2, x3) - calculate_G(x1 - h, x2, x3)) / (2 * h)
    dG2 = (calculate_G(x1, x2 + h, x3) - calculate_G(x1, x2 - h, x3)) / (2 * h)
    dG3 = (calculate_G(x1, x2, x3 + h) - calculate_G(x1, x2, x3 - h)) / (2 * h)
    # # 验证

    # Calculate deltaGB values
    deltaGB1 = G + dG1 - x1 * dG1 - x2 * dG2 - x3 * dG3
    deltaGB2 = G + dG2 - x1 * dG1 - x2 * dG2 - x3 * dG3
    deltaGB3 = G + dG3 - x1 * dG1 - x2 * dG2 - x3 * dG3

    # print(deltaGB1, deltaGB2, deltaGB3)

    # 北科博士段论文写法，去除第一组元偏导的项
    GB1 = G -x2*dG2 -x3*dG3
    GB2 = G +(1-x2)*dG2 -x3*dG3
    GB3 = G -x2*dG2 +(1-x3)*dG3

    # print(GB1, GB2, GB3)

    # 以上是计算超额吉布斯能的，以下通过多元butler方程计算表面张力

    R = 8.314
    N0 = 6.02E+23
    # 二元 体系
    L = 1.091
    A = [L * pow(N0,1/3) * pow(M[0]/rou[0], 2/3), L * pow(N0, 1/3) * pow(M[1]/rou[1], 2/3), L * pow(N0, 1/3) * pow(M[2]/rou[2], 2/3)]

    def equations(vars):
        xs1, xs2 = vars
        if xs1 <= 0 or xs1 >= 1 or xs2 <= 0 or xs2 >= 1 or (1 - xs1 - xs2) <= 0:
            return [1e6, 1e6]  # 返回大数以表示无效解
        try:
            term1 = theta[0] + (R * T * log(xs1 / x1) + (beta - 1) * GB1) / A[0]
            term2 = theta[1] + (R * T * log(xs2 / x2) + (beta - 1) * GB2) / A[1]
            term3 = theta[2] + (R * T * log((1 - xs1 - xs2) / x3) + (beta - 1) * GB3) / A[2]

            # 方程组
            eq1 = term1 - term3
            eq2 = term2 - term3
            return [eq1, eq2]
        except (ValueError, ZeroDivisionError):
            return [1e6, 1e6]  # 返回大数以表示无效解

    # initial_guess = [0.2, 0.2]
    # solution = fsolve(equations, initial_guess)
    # xs1_solution, xs2_solution = solution
    #
    # solution = root(equations, initial_guess, method='lm')  # 使用Levenberg-Marquardt算法
    from scipy.optimize import root
    initial_guess = [0.1, 0.1]
    solution = root(equations, initial_guess, method='lm')  # 使用Levenberg-Marquardt算法

    if solution.success:
        xs1_solution, xs2_solution = solution.x
        result = theta[0] + (R * T * log(xs1_solution / x1) + (beta - 1) * GB1) / A[0]
    else:
        print("优化未收敛，请检查初始条件")

    return theta[0] + (R * T * log(xs1_solution / x1) + (beta - 1) * GB1) / A[0]

def calculate_binary_viscosity(system,T,x1,x2):

    L12,V = parameter_vicosity(system,T)

    V_PE = x1*x2*(L12[0]+L12[1]*(x1-x2)+L12[2]*(x1-x2)**2+L12[3]*(x1-x2)**3)

    return x1*V[0]+x2*V[1]+V_PE

def calculate_ternary_viscosity(system,T,x1,x2,x3):

    elements = system.split('_')
    system12 = elements[0]+'_'+elements[1]
    system23 = elements[1]+'_'+elements[2]
    system13 = elements[0]+'_'+elements[2]

    L12,V12 = parameter_vicosity(system12,T)
    L23,V23 = parameter_vicosity(system23,T)
    L13,V13 = parameter_vicosity(system13,T)

    V = [V12[0],V12[1],V23[1]]

    V_i = x1*V[0]+x2*V[1]+x3*V[2]

    calculator = PE_cal.PE_calculator()
    V_PE = calculator.GSM(x1,x2,x3,L12[0],L12[1],L12[2],L12[3],L13[0],L13[1],L13[2],L13[3],L23[0],L23[1],L23[2],L23[3])

    return V_i+V_PE